DTI Crossing Fibers Estimation Using Fast 
Independent Component Analysis* 

Ch. Khodaverdian, H. Soltanian-Zadeh Senior Member, IEEE 



Abstract — Tracing fibers in voxels containing multiple fiber 
tracts has been an important problem in white matter fiber 
tractography using Diffusion Tensor Magnetic Resonance 
Imaging (DT-MRI) data. In order to solve this problem, based 
on the fact that DTI signals have super-Gaussian pdf's, fast 
Independent Component Analysis (ICA) is applied to 
decompose the original DTI signal from voxels containing 
multiple fibers into its non-Gaussian components. 
Tractography is then done based on the eigenvectors of the 
new components. In this approach, the voxels with crossing 
fibers are identified based on the orientation difference, 
number of fibers coexisting in a voxel, and the q value. 
Results show more reconstructed fibers by the proposed 
method compared to the previous methods. 



I. Introduction 



D 



T t 



(D II) 



leioiuct In aging (M II) technique till provides 
inform alioi about tie random i otioi of w iter i oltcnles 
referred to is B row liji % otioi. Tit w iter diffusion is 
anisotropic ii lit hair, i lilt n itter. In fact, w iter tends to 
diffuse predoi iutely ii i dirtttioi parallel to tie axons, 
so His anisotropy en reveal i icroscopic properties of tie 
niton j of He nerve filers. I lerefore, I) I I t a i 1 e i se i to 
no-invasivtly assess i lilt i alter fiber pathi ajs aid also 
to build connectivity maps [1], 

One of tie major concerns in in -vivo w lite i alter 
fiber tractography is He tracing of tic fibers in voxels 
containing i iltiplt fibtrs w lib different o rit n la tio n s. In 
conventional stream line tracking m etbods, lie direction of 
He eigenvector corresponding to He largest eigenvalue of 
tie diffusion tensor is tonsidtred to be He local direction 



Excellence. Electrical and Computer Engineering Departmenl. Universih 
of Tehran, Tehran 14395-5 15. Iran (email: christin mkh@ vahoo.com ). 

H. Soltanian-Zadeh is with the Control and Processing Center of 

II 11 Ekclncal aikl il)| ler Enjin nil D rlillcill. I nivt 1 it; 

of Tehran, Tehran 14395-515, Iran and Image Analysis Laboratory, 

Radiology Department Henry Ford Health System, Detroit MI, 48202 

USA (emails: hszadeh@ut.ac.ir ; hamids@rad.hfh.edu ). 

"This work was supported in part by the Iran Telecommunication 
Research Center (ITRC), Tehran, Iran. 



of the fiber pathways [l]-[2]. However, this direction is not 
valid in voxels with multiple fibers. 

In order to solve this problem, we should determine which 
voxels contain crossing voxels. It is known that in voxels 
containing non-parallel fibers, the linear anisotropy ( C, ) is 

low but parallel anisotropy ( C ) is high. However, 
determination of an appropriate threshold for C is 
difficult. In addition, the position where fibers cross affects 
the C p value [3]. 

As can be seen in Fig. 1, a crossing has occurred but the 
C value is low. So, using the C value will not be 
sufficient for accurate classification of crossing voxels. 

In this paper, we consider a voxel to be a crossing 
voxel if: 1) multiple fibers coexist in that voxel; and 2) the 
C p value is sufficiently high. 

We employ conventional Independent Component 
Analysis (ICA) to decompose multiple tensors for the 
voxels identified as containing crossing fibers. This method 
is called "fast ICA" and finds multiple non-Gaussian 
sources from multi-channel data by maximizing their non- 
Gaussianity [4]. 



II. METHODS AND MATERIALS 

A. Determining Crossing \'oxels 
We identify a voxel with crossing fibers based on the 
following criteria: 1) multiple fibers coexist in the voxel; 
and 2) the orientation difference in the previous voxel 
exceeds a certain degree [3]; and C is larger than a 
threshold. According to the fact that white matter fibers 
have low curvatures, two fibers with large orientation 
differences in the previous voxel maintain their difference 
when they meet in the so called crossing voxels. 

To this end, we ran the streamline tracking algorithm in 
the whole brain. Then, a Region of Interest (ROI) was 
selected and each of the voxels within this (ROI) was 
checked whether they were contained in multiple fibers or 
not. If yes, the fibers containing the selected voxels were 
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listed and the orientation difference between these fibers in 
a step before reaching the desired voxel was calculated. 

When the maximum orientation difference between all 
possible pairs of fibers passing through a voxel exceeds a 
specified threshold, this voxel will be defined as a crossing 

The orientation difference of the two fibers is calculated 
as [3]: 
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where V i corresponds to the i-th eigenvector. 
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Fig. 1 . Position of crossing fibers affects the value of C p . For (b) w. 



expecl a high C value i 



slowL3j. 



= log 



s(gj,b) 
s(gj,b = 0) 



(2) 





% 




hi 


hi 


hi 




ha 




hi 


^ 


hi 


h; 


h«i 


hu 


*iy 


hij 


% 




hi 




hii 


hij 


hM 




hij 





Fig. 2. The cubic window used for the 

in 3-D neighborhood of the centra] voxel |4J. 

After construction of the measured data matrix X, fast 
ICA was used to extract a new set of statistically 
independent vectors of Y called independent components 
[6]: 



Y=WX 



(3) 



Among the selected voxels, those having C value less 
than a selected threshold were discarded. In previous 
studies, the threshold was set at 0.3 [4], while in this work, 
the threshold is set at a lower level. The reasons for the 
above setting were consideration of the pervious two 
criteria and selection of crossing voxels with low 
C value. 

Having identified the crossing voxels, it is time to 
decompose their DTI signals into non-Guassian 
components. 

B. Decomposition of DTI Signal 
According to the fact that the width of a fiber tract may be 
more than the dimension of a voxel, crossing fibers can 
cover several neighboring voxels [4]-[5]. In the proposed 
approach, it is assumed that it covers 19 adjacent voxels 
(see Fig. 2) [4]. 

The data matrix X was constructed for the cubic window 
described above. Each element of X, x v , correspond to i-th 
voxel and j-th gradient DTI signal and is defined by the 
logarithm of the division of the j-th gradient signal 
s{g ,b) by the j-th gradient signal s(gj,b = 0) [4], 



where W is the unmixing matrix whose columns contain 
unmixing weights corresponding to individual sources [6]. 
Here, it was assumed that the source signals were 
decomposed into two independent components using a 
deflationary orthogonalization method that finds the 
components sequentially. Each of the components 
corresponds to one of the crossing fibers. 



C. Fiber Tracking 
The two estimated independent components were first 
projected onto the data space. The first and second 
components of the central voxel (10-th component) were 
used to construct two diffusion tensors. In the next step, the 
eigenvectors and eigenvalues of the constructed diffusion 
tensors were calculated and then fiber tracking was 
performed using FACT (Fiber Assignment by Continuous 
tracking) algorithm. 

In the FACT algorithm, the 3-D trajectories are 
reconstructed from a 3-D vector field by propagating a line 
from a seed point by following the local vector orientation. 
The propagation direction is alerted at voxel boundary 
interfaces. Line propagation is terminated according to the 
following 2 criteria: 1) The extent of anisotropy. In low 
anisotropy regions, such as gray matter, the fractional 
anisotropy is in the range of 0.1-0.2. In such regions the 
fastest diffusivity direction is not well defined and the 
largest principle axis is prone to noise error. So the FA 
threshold for propagation termination can be set to 0.2. 2) 
Large angular change. In diffusion tensor calculation, it is 
assumed that there is no sharp turn during line propagation. 
We have considered the angle-change threshold to be 
45° [1]. 
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The tractography was done twice, once selecting the ICA Table I illustrates the quantitative difference between the 
component corresponding to the larger primary eigenvalue two methods, 
and once for the other component. 

III. EVALUATION OF FAST ICA 

In order to evaluate the performance of fast ICA, DT-MRI 
data were simulated for two crossing fibers which cover 
five voxels with different mixing ratios in each voxel. The 
angle between the crossing fibers were considered to be 45 
degrees. 

The estimated fiber tracts can be seen in Fig. 3. The 
mixing ratio of the two horizontal (black bar) and diagonal 
(gray bar) tracts, considered in this estimation were 
0.5:0.58 in the central voxel and 0:0.58, 0.5:0.06, 0.5:0.06 
0:0.58, in four adjacent voxels [4]. ICA was applied to 
data matrix X constructed based on the estimated fiber 
tracts. The mean error in 100 trails was 12.8° which 
showed almost reliable performance of fast ICA. 




IV. DATA ACQUISITION 

The DTI data used in this experiment were acquired using a 
1.5 Tesla, 63.8859 MHz, GE Signa System (GE Medical 
Systems, Milwaukee, WI), 256x256 matrix, 29, 5mm thick 
slices, 25 gradient directions with b=1000sl mm 2 and 
one b = acquisition. 

V. EXPERIMENTAL RESULTS 

Results of the ICA decomposition are shown in Fig. 4. The 
selected ROI is indicated by a blue box in Fig. 4(a). This 
region is likely to contain several crossing fibers. The blue 
arrows in Fig. 4(b) correspond to the primary eigenvectors 
of voxels in the selected ROI. The red and green arrows 
represent the primary eigenvectors of the two independent 
components corresponding to the voxels identified as 
crossing voxels, with maximum angular difference > 20 
degrees and C value > 0.08. 

The fiber tractography results for the selected ROI are 
shown in Fig. 5. Fiber tracts represented in Fig. 5(a) and 
Fig. 5(b) are obtained using conventional FACT tracking 
algorithm where Fig. 5(c) and Fig. 5(d) are the 
reconstructed fibers using FACT algorithm by considering 
two independent components. As seen in the figures, the 
second method results in tracking of many more fibers. 
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VI. CONCLUSIONS 

White matter fiber tractography in voxels containing 
multiple fibers is a difficult problem. In solving the 
problem, the main concerns are: 1) How to identify the 
voxels containing the crossing fibers; 2) How to perform 
the fiber tractography. 

In this work, for the identification of the crossing voxels, 
in addition to the C value, coexistence of multiple fibers 
and their orientation differences were considered as well. 
Then, fast ICA method was used to decompose the DTI 
source signals of the selected voxels into two non-Gaussian 
components. Next, a FACT tractography algorithm was 
preformed twice, once for each of the independent 
components. Although experimental results showed better 
reconstruction of fiber tracts, cross-validating the results 
are not yet possible via noninvasive methods. 

In this study, we limited the number of crossing fibers to 
two and estimated two independent components each 
corresponding to one of the crossing fibers. However, since 
ICA works as long as the number of independent sources to 
be estimated is less than the number of input channels, and 
we have nineteen input channels in this study, additional 
fibers within a voxel can be detected in the future work. 
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